1. Field of the Invention
The invention relates to a method for reconstructing a three-dimensional image of an object or a subject scanned, preferably in linear or circular fashion, in the context of a tomosyntheisis procedure of the type wherein during the scanning, several individual projection images located in a 2D projection image space are recorded in the form of digital projection image data of the object located in a 3D object space, and are back-projected into a 3D reconstruction image volume in order to produce the reconstruction image, the object being irradiated with X-rays from various projection angles xcfx86 for the purpose of recording the projection images, and wherein the radiation exiting from the object is recorded using a detector that supplies digital output image signals, the output image signals representing the projection image data and being supplied to a computer for image reconstruction.
2. Description of the Prior Art
In medical imaging systems, it is known to record an object tomosynthetically and to reconstruct it three-dimensionally. Since in the tomosynthetic data recording the object to be imaged are projected onto a detector from only a few spatial directions, the object scanning is incomplete. This is expressed in a poor depth resolution of the 3D reconstruction. In a simple back-projection (identical to the summation of displaced projection images, classical slice method), without reconstructive corrections locus-frequency-dependent artefacts are contained in the tomograms. Given larger objects, these disturbances can be very extensive.
The following demands are placed on a good 3D reconstruction: best possible suppression of structures foreign to the slice, defined slice behavior or characteristic, i.e. depth resolution independent of locus frequency, purposive controlling of the characteristics of the reconstructed tomogram.
In system-theoretical terms, the tomosynthetic imaging of an object distribution f(x, y, z) to form a 3D reconstruction image g(x, y, z) can be formulated as a convolution (indicated, as is standard, by the symbol *) of the object distribution with the point image function h(x, y, z) of the imaging process:
g(x, y, z)=h(x, y, z)*f(x, y, z)
The point image function h(x, y, z) describes both the xe2x80x9cmeasurement processxe2x80x9d (projection and back-projection) and also reconstructive measures, such as e.g. filterings. In the 3D Fourier space, the Fourier-transformed image G(xcfx89x,xcfx89y,xcfx89z) is described by a multiplication of the Fourier-transformed object distribution F(xcfx89x,xcfx89y,xcfx89z) with the 3D transmission function (or Modulation Transfer Function MTF) H(xcfx89x,xcfx89y,xcfx89z):
G(xcfx89x,xcfx89y,xcfx89z)=H(xcfx89x,xcfx89y,xcfx89z)xc2x7F(xcfx89x,xcfx89y,xcfx89z)
with F(xcfx89x,xcfx89y,xcfx89z) as the 3D Fourier transformation of the locus distribution f(x, y, z)
F(xcfx89x,xcfx89y,xcfx89z)=FzFyFxf(x, y, z)
and analogously for the image g(x, y, z) and the point image function h(x, y, z). By modification of the point image function h(x, y, z) or of the Modulation Transfer Function H(xcfx89x,xcfx89y,xcfx89z) the imaging process can be purposively influenced.
The image quality of a reconstructed tomogram can be judged using the slice transmission function h(xcfx89x,xcfx89y|z) [D. G. Grant, TOMOSYNTHESIS: A Three-Dimensional Radiographic Imaging Technique, IEEE Trans. on Biomed. Eng. 19 (1972), 20-28].
h(xcfx89x,xcfx89y|z)=Fzxe2x88x921 H(xcfx89x,xcfx89y,xcfx89z)
The slice transmission function h(xcfx89x,xcfx89y|z) is a hybrid of a representation in the Fourier space and in the locus space. It indicates the spectral content (frequency response) with which objects at a distance z from the reconstruction plane contribute to the reconstruction image. For z=0, it indicates the spectral content of the slice to be reconstructed, and for small z it indicates the screen characteristic of a reconstruction slice of finite thickness, and for large z, specifically in tomosynthesis, it indicates the locus-frequency-dependent feed-through of artefacts. The desired image characteristics (see above) can be formulated as follows using the slice transmission function.
A good suppression of structures foreign to the slice signifies a rapid decay of the slice transmission function h(xcfx89x,xcfx89y|z) into z.
Defined slice behavior, i.e. a slice profile independent of the locus frequency, is achieved by a separation of the slice transmission function h(xcfx89x,xcfx89y|z) into a portion Hspectrum(xcfx89x,xcfx89y), which represents the spectral image content of the reconstructed tomogram, and into a portion hprofile(z) that defines the slice profile. A complete separation is in general not possible in tomosynthesis, since as the locus frequencies decrease the scanning becomes increasingly incomplete. It is already a considerable advantage, however, if a separation is achieved for a limited locus frequency region:
xe2x80x83h(xcfx89x,xcfx89y|z)=Hspectrum(xcfx89x,xcfx89y)xc2x7hprofile(z)
In the literature, essentially the following methods of reconstruction in tomosynthesis are found:
Simple back-projection: [D. G. Grant, TOMOSYNTHESIS: A Three-Dimensional Radiographic Imaging Technique, IEEE Trans. on Biomed. Eng. 19 (1972), 20-28]: As in the classical slice method, by simple summation of the projection images an uncorrected reconstruction image is obtained. The method is rapid, simple and robust, but yields poor image results. Simple heuristic 2D filterings of the reconstructed tomograms [R. A. J. Groenhuis, R. L. Webber and U. E. Ruttimann, Computerized Tomosynthesis of Dental Tissues, Oral Surg. 56 (1983), 206-214] improve the image impression, but do not provide image material that can be reliably interpreted, and are thus not satisfactory.
Influencing of the recording geometry: By the selection of the sampling curve, such as e.g. spirals or threefold concentric circular scanning [U. E. Ruttimann, X. Qi and R. L. Webber, An Optimal Synthetic Aperture for Circular Tomosynthesis, Med. Phys. 16 (1989), 398-405], the point image function can be manipulated. The methods require a high mechanical outlay and are inflexible, since the image corrections are already determined with the measurement data. In addition, only slight image improvements have been achieved.
Iterative methods: [U. E. Ruttimann, R. A. J. Groenhuis and R. L. Webber, Restoration of Digital Multiplane Tomosynthesis by a Constrained Iteration Method, IEEE Trans. on Medical Imaging 3 (1984), 141-148]. With the aid of an iterative reconstruction, in principle the measurement process can be flexibly modeled and corrected by approximation. The iterative method attempts to invert the point image function. Since this is not possible, due to the incomplete scanning, secondary conditions must be introduced in order to guarantee the unambiguity of the solution. The methods achieve good image quality, but are not obvious and are difficult to use. In particular, the required introduction of secondary conditions can strongly falsify the image impression in an undesired manner. In addition, there are problems with the stability of the algorithms, and the computing time is prohibitive for a routine application.
Algebraic methods: For a set of tomograms, the point image function is set up slice-by-slice, which leads to a matrix formulation of the point image function. The exact inversion of the matrix is not possible, due to the incomplete scanning in the tomosynthesis. [J. T. Dobbins, Matrix Inversion Tomosynthesis Improvements in Longitudinal X-Ray Slice Imaging, U.S. Pat. No. 4,903,204, Feb. 20, 1990] falsely claims invertibility. This method, however, uses multiple subsequent corrections of the reconstruction images, which underlines the inadequacy of his approach. No publication of the image results, however, is available.
Reconstructive filtering methods: The present application is based on this principle. In the literature various approaches can be found, which can be classified on the basis of the space in which the filtering takes place. From the theoretical point of view, filtering in the 3D space of the object imaged or in the 2D space of the projection images are equivalent. The filters to be used can be transformed correspondingly.
Filtering in the 3D Fourier space of the imaged object [H. Matsuo, A. Iwata, I. Horiba and N. Suzumura, Three-Dimensional Image Reconstruction by Digital Tomosynthesis Using Inverse Filtering, IEEE Trans. on Medical Imaging 12 (1993), 307-313]: A filtering function is defined in the 3D Fourier space. A 3D reconstruction volume is obtained as a set of tomograms by simple back-projection. The reconstruction volume is transformed into the locus frequency space with the aid of the rapid Fourier transformation, and is filtered in that space and is back-transformed again. The method is essentially more computing-intensive than the filtering at the projection images, and is sensitive to discretization and to the size of the filtered reconstruction volume. From the theoretical point of view, [H. Matsuo, A. Iwata, I. Horiba and N. Suzumara, Three-Dimensional Image Reconstruction by Digital Tomosynthesis Using Inverse Filtering, IEEE Trans. on Medical Imaging 12 (1993), 307-313] have indeed calculated the correct 3D transmission function of the circular tomosynthesis function. The modification of the inverse filtering function which is selected mirrors the relations of a complete object scanning (computed tomography (CT) with surface detector and parallel beam approximation). The filtering function used corresponds to the generally known Shepp-Logan CT filter. The approach borrowed from CT does not do justice to the relations of tomosynthesis, and the image quality cannot be optimal.
Filtering in the 2D Fourier space of the projections (filtered back-projection): In the 2D Fourier space of the projections, a filtering is defined. Using fast Fourier transformation, the projection images are successively transformed into the locus frequency space, are filtered in that space, are back-transformed, and are subsequently calculated to form a reconstruction image by means of back-projection. The filtering can also be carried out in equivalent fashion in the locus region, by means of 2D convolution. The method is essentially faster than the 3D filtering, and has no artificial instabilities due to method parameters. The present invention likewise belongs to the class of filtered back-projections. From the literature, only ectomography [P. Edholm, G. Graniund, H. Knutsson and C. Petersson, Ectomographyxe2x80x94A New Radiographic Method for Reproducing a Selected Slice by Varying Thickness, Acta Radiologica 21 (1980), 433-442] is known as a representative. In ectomography, however, the 2D filtering function is set up only empirically. The filters cannot be designed purposely for a particular optimization task, and an image manipulation, and thus a high image quality, is thus not achieved.
An object of the present invention is to provide a method for image reconstruction in which the image quality of the reconstruction image can be manipulated, and is improved in relation to the prior art.
The present invention provides the use of 2D filters in the filtered back-projection of the tomosynthesis for the production of the reconstruction image. On the basis of an analysis of the 3D transmission function, a filtering function is produced and is applied in the computing means for the reconstruction of the image. The entire procedure is divided into seven steps:
1. Calculation of a 3D transmission function Hproj(xcfx89x,xcfx89y,xcfx89z) from the recording geometry for the individual projection image recording, and the back-projection of the individual projection images into the 3D reconstruction image volume,
2. Approximate inversion of the 3D transmission function Hproj(xcfx89x,xcfx89y,xcfx89z) for the determination of an inversion function Hinv(xcfx89x,xcfx89y,xcfx89z),
3. Production of a 3D filtering function Hopt(xcfx89x,xcfx89y,xcfx89z) dependent on one or more desired image characteristics of the reconstruction image,
4. Determination of a resulting 3D filtering function Hfilter(xcfx89x,xcfx89y,xcfx89z) by multiplication of the 3D filtering function Hopt(xcfx89x,xcfx89y,xcfx89z) and the inversion function Hinv(xcfx89x,xcfx89y,xcfx89z),
5. Determination of a 2D filtering function H2Dfilter,xcfx86(xcfx89u,xcfx89v) from the resulting 3D filtering function Hfilter(xcfx89x,xcfx89y,xcfx89z), by coordinate transformation of the 3D object space into the 2D projection image space of the respective individual projection images, at the projection angle xcfx86,
after which the reconstruction of the image takes place in the computing means, with the following steps:
6. Application of the 2D filtering function H2Dfilter,xcfx86(xcfx89u,xcfx89v) to the associated individual projection image data,
7. Production of the reconstruction image by back-projection of the individual projection image data, filtered according to 6., into the 3D reconstruction image volume.
The overall 3D transmission function thus results from the product of the individual components:       H    ⁡          (                        ω          x                ,                  ω          y                ,                  ω          z                    )        =                                          H            opt                    ⁡                      (                                          ω                x                            ,                              ω                y                            ,                              ω                z                                      )                          ·                              H            inv                    ⁡                      (                                          ω                x                            ,                              ω                y                            ,                              ω                z                                      )                                      ⏟                  3          ⁢          D          ⁢                      xe2x80x83                    ⁢          filtering          ⁢                                    xe2x80x83                        ⁢                          xe2x80x83                                ⁢          function          ⁢                      xe2x80x83                    ⁢          to          ⁢                      xe2x80x83                    ⁢          be          ⁢                      xe2x80x83                    ⁢          applied          ⁢                      xe2x80x83                    ⁢                                    H              filter                        ⁡                          (                                                ω                  x                                ,                                  ω                  y                                ,                                  ω                  z                                            )                                            ·                  H        proj            ⁡              (                              ω            x                    ,                      ω            y                    ,                      ω            z                          )            
The object of the invention is also achieved in a procedure divided into six steps:
1. calculation of a 3D transmission function Hproj(xcfx89x,xcfx89y,xcfx89z) from the recording geometry for the individual projection image recordings, and back-projection of the individual projection images into the 3D reconstruction image volume,
2. Approximate inversion of the 3D transmission function Hproj(xcfx89x,xcfx89y,xcfx89z) for the determination of an inversion function Hinv(xcfx89x,xcfx89y,xcfx89z),
3. Determination of a 3D filtering function Hopt(xcfx89x,xcfx89y,xcfx89z) dependent on one or more desired image characteristics of the reconstruction image,
4. Determination of a resulting 3D filtering function Hfilter(xcfx89x,xcfx89y,xcfx89z) by multiplication of the 3D filtering function Hopt(xcfx89x,xcfx89y,xcfx89z) and the inversion function Hinv(xcfx89x,xcfx89y,xcfx89z),
after which the reconstruction of the image takes place in the computing means with the following steps:
5. Production of the reconstruction image by back-projection of the individual projection image data into the 3D reconstruction image volume,
6. Application of the resulting 3D filtering function Hfilter(xcfx89x,xcfx89y,xcfx89z) to the reconstruction image for the filtering thereof.
Here, the reconstruction image itself is filtered, after its production, with the resulting 3D filtering function, in contrast to the solution previously described, in which the individual projection images are already filtered before the back-projection. Both methods have in common the application of the 3D filtering function according to the desired image characteristics, so that the reconstruction image shows the desired image characteristics.
The invention is suitable for all radiological applications, such as imagings of the skull, of the skeletal system and of the inner organs, as well as for applications in nuclear medicine with specific collimators.